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Simulations of two-planet systems through all phases of stellar 
evolution: implications for the instability boundary and white 
dwarf pollution 
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ABSTRACT 

Exoplanets have been observed at many stages of their host star's hfe, including the main se- 
quence (MS), subgiant and red giant branch stages. Also, polluted white dwarfs (WDs) likely 
represent dynamically active systems at late times. Here, we perform 3-body simulations 
which include realistic post-MS stellar mass loss and span the entire lifetime of exosystems 
with two massive planets, from the endpoint of formation to several Gyr into the WD phase 
of the host star. We find that both MS and WD systems experience ejections and star-planet 
collisions (Lagrange instability) even if the planet-planet separation well-exceeds the analyti- 
cal orbit-crossing (Hill instability) boundary. Consequently, MS-stable planets do not need to 
be closely-packed to experience instability during the WD phase. This instability may pollute 
the WD directly through collisions, or, more likely, indirectly through increased scattering of 
smaller bodies such as asteroids or comets. Our simulations show that this instability occurs 
predominately between tens of Myr to a few Gyrs of WD cooling. 

Key words: planet-star interactions, planets and satellites: dynamical evolution and sta- 
bility, stars: evolution, stars: AGB and post-AGB, stars: white dwarfs 



1 INTRODUCTION 

A planet's life may be split into four distinct stages: 1) formation 
and concurrent dynamical excitation, 2) main sequence (MS) 
evolution, 3) evolution during post-MS stellar phase changes, 
and 4) white dwarf (WD) evolution. The first stage generally 
lasts no longer than 0.1% of the entire MS lifetime. The second 
stage is relatively dynamically quiescent, with only occasional 
but often important scattering interactions. In the third stage, 
the planet is subject to dynamical changes due to the star's vi- 
olent actions as it becomes a giant. In the final stage, the star 
has become a WD, and the planet again enters and remains in a 
phase of relative dynamical quiescence occasionally punctuated 
by scattering interactions or external forcing. This general pic- 
ture, which does not include possibilities such as the capture of 
free-floating planets, planetary destruction due to supernovae, 
or multiple host stars, describes the life cycle of the vast major- 
ity of known exoplanets. 

The volume of planetary literature investigating the 



* E-mail: veras@ast.cam.ac.uk 

t E-mail: alex.mustill@uam.es 

X E-mail: amy.bonsorOobs.ujf-grenoble.fr 

§ E-mail: wyatt@ast.cam.ac.uk 



first two stages dwarfs the literature describing the fi- 
nal two stages, despite t he fact that the Universe is al- 
ready over 13.5 Gyr old (|jarosik et al.l |2011^ and that the 
Milky Way contains about 10^ WDs l|Binnev fc Tremaind 
l2008l . Pgs. 2-3 and iHolberg et al] |2008| ). Further, these fi- 
nal two stages are becoming increasingly relevant given the 
sugge stions or d iscoveries of exoplanet s in post-MS systems 



Lee et al.l [20091: iMullallv et al.| l2009l: ISetiawan et al 



Wickramasinghe et al.l 2010l: 



Charpinet et al.l 



2011 



Adamow ct al.' ' 2OI2I : iFarihi et all l2012al : iLee et all l2012al ll:: 
Sat o ct al. 2012^). 

Explorations of exosyste m evolution in the thir d stage 



include one- plan e t studies (IVifla,ver fc Livid |2007|. |2009|: 



2012 



Veras et al.l I2OIII: iKratter fc PeretsI 120 121: iMustill fc Viflaverl 



Nordhaus fc Spiegelll2012l: ISpiegel fc Madhusud haDll2012l : 



Veras fc Tout! I2OI2I : lAdams et al.l |2013|). just 



few ded- 
icated multiple-plane t studies llDebes fc SigurdssorJ |2002| : 
iPortegies I2OI2I : IVovatzis et all |2013|), and studies fo- 
cusing on 



the evolut ion of comets jAlcock et al. Il986l : 



iParriott fc AlcocklflQOl ) . Further. iBonsor fc Wvatd l|2010l ') con- 
sider the effect of post-MS evolution on debris discs. Motivated 
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by o bserv ations of metal-pol luted WDs, I Bonsor et al ] (|201ll . 
l2012f ) and iDebes et alj (|2012l ') model the interplay between a 
planet and a belt of smaller material amidst stellar mass loss. 

Here, we self-consistently simulate the second, third and 
fourth stages together. We combine stellar evolution with plan- 
etary gravitational scattering amongst multiple massive planets, 
and extend the work of Dcbes & Sigurdsson (2002) by consider- 
ing full-lifetime simulations with realistic mass loss prescriptions 
at each post-MS phase. Following the evolution over the whole 
stellar lifetime means that we can be sure that systems whose 
stability is investigated on the giant and WD stages will have 
survived the long MS evolution. Through these integrations, we 
can determine what types of planetary architectures might be 
expected in exoplanet-hosting WD systems, and could allow us 
to extrapolate backwards in time from observed WD systems. 
We restrict our explorations to two-planet systems in this initial 
study given the vast phase space to explore; three-planet simu- 
lations will be presented in a follow-up paper. First, we briefly 
summarize our knowledge of planetary instability for one- and 
two-planet systems during the MS (Subsection 1.1) and post-MS 
(Subsection 1.2). 



colliding with the star. If both planets remain bound and re- 
tain their ordering, and no collision with the star occurs, then 
the system is "Lagrange stable'Q- Unlike Hill stability, Lagrange 
stability does not benefit from a known analytical formulation, 
but rather empirical estimates based on numerical simulations. 

The analytical Hill stability boundary is conservative. Two 
planets whose initial separation is less than the Hill stable dis- 
tance may in fact remain stable. If the initial separation is 
greater than the Hill stable distance, then the planets are guar- 
anteed to retain their ordering. In simulations of th e HD 12661 
and 47 Uma systems, iBarnes fc Greenber j |2006^ found that 
pairs of planets close to the Hill stability boundary are not La- 
grange stable, and hence are not generally stable. They tenta- 
tively suggest that the Lagrange stability boundary exceeds the 
Hill stable boundary by at least 21 % as measured by the semi - 
major axis ratio. Subsequent work (Barnes fc Greenberg||2007l ) 
revealed how mean motion commensurabilities can broaden 
the divide between t h e Hil l and Lagrange stable boundaries; 
iKopparapu fc Barnes! (|201G| ) demonstrated how the boundary 
between stable and unstable systems is not sharp. Hence, nu- 
merical validation of analytical stability estimates is crucial. 



1.1 Instability in Main Sequence Planetary Systems 

Dynamical instability in planetary systems is often said to occur 
when a planet suffers a close encounter with the star or another 
planet, or is ejected from the system. Occasionally, investigators 
use stricter definitions of instability, such as when the semimajor 
axis or eccentricity variation of a planet exceeds a certain per 
cent of its nominal value. Additionally, a wide body of literature 
has ar isen characteri z ing ch aotic orbits as a precursor to insta- 
bility; iDarriba et al.l l|2012l ) and references therein summarize 
many of these techniques. 



1.1.1 One-Planet Instability 

One planet orbiting a MS star will typically remain stable 
throughout the star's MS lifetime in the absence of external 
forces. Exceptions may include planets which a re close enough 



to their parent stars to be tid ally disrupted (e.g. | Gu et al.|[2003l ) 
and possibly evaporated (e.g. iGuillot et ai]|l996h . In the oppo- 
site extreme, a planet which is far enough away from its parent 
star may be ejected due to external forces such as passing stars 
(Zakamska fc Trcmainc 2004; Vcras & Mocckcl 2012) or achieve 
a high enough eccen tricity through Galacti c tides to cause a col- 
lision with the star (|Veras fc Evans|[2013al lbh. 



1.2 Instability in Post-Main Sequence Planetary 
Systems 

Mass loss from a dying star can trigger planetary instability in 
different ways, which are outlined below. A common assump- 
tion amongst the studies which have considered instabilities in 
post-MS systems is isotropic stellar mass loss. We also adopt this 
assumption here, as modelling non-isotropic mass loss would sig- 
nificantly complicate both numerical and analytical modelling 
and is best left to sep arate, dedicat ed studies. One such dedi- 
cated post-MS study (Parriott fc Alc ock 1998) importantly ob- 
serves that the speed of (effectively massless) comets near the 
boundary of a planetary system may be comparable to the re- 
coil velocity of the parent star due to asymmetric mass loss. 
That study suggests that anisotropic mass loss will affect the 
details of planets being ejected after scattering but is unlikely 
to have a significant effect on the prior dynamics. Other studies 
modelling planetary dynamics due to non-isotropic mass loss 
instead focus on jet ac c elerations p resent at the birth sites of 
planets (|Namouniil2005l . |2007| . [20il ). 

Further, in all cases we assume the planets are or- 
biting a single star. Extensions to the mu lt iple-star case 
(IKratter fc Peretd |2012| : iPortegies ZwartI |2012| : IVeras fc Tou^ 
I2OI2I ) are likely to be nontrivial. 



1.1.2 Two-Planet Instability 

In addition to tidal interactions and external forces, the mutual 
perturbations between two planets may also create instability. 
Partially motivated by tractable analytical solutions to the gen- 
eral three-body problem, the source of this instability has been 
studied extensively. If the orbits of two planets are guaranteed 
to never overlap (precluding a collision between bo t h pla nets), 
then they are said to be "Hill stable". 'Gladman' (imi) pio- 
neered the analytic use of Hill stability for planetary systems in 
specific cases and has motivat ed many subsequent an al yses, a s 
recently summarized by, e.g.. iDonnisoiil (|2009l , l2010al lbl 120111 ). 
Hill stability does not guarantee that the outer planet remains 
bound to the system, nor does it prevent the inner planet from 



1.2.1 One-Planet Instability 

For decades, binary star investigations revealed that stellar mass 
loss causes orbital semimajor axis expansion. Less well-known 
is that when the mass loss is rapid enough, t he eccentricity 
of the companion's orbit can change as weU (|Omarovl [l96a : 
iHadiide mctriou 1963, 1966; Vcras ot al. 201l[). If the eccentric- 
ity is great enough, a planetary companion may escape from the 
system. 

^ T his type of stability has also been referred to as "Laplace stability" 
(e.g. lKubala et al.ll993h and featured but remained unnamed in many 
papers published before the discovery of exoplanets and high-speed 
computing. 
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In the more commonly-used adiabatic limit, de/dt = and 
dajdt = - (flip,) {dp/dt). Here, p = G (M* + Mp), and Mp 
are the masses of the star and planet, a is the planet's semimajor 
axis, and e is the planet's eccentricity. This limit holds when the 
mass loss timescale is much longer than the orbital timescale. 
Adiabaticity will be broken, even briefly, if at any point a sudden 
burst of mass loss causes the timescales to become comparable 
l|Veras fc Toutll2012l : IVeras fc Wvattll2012l ). Hence, characteriz- 
ing whether or not planetary evolution is adiabatic amidst mass 
loss will be important for any post-MS scattering study. 

Another source of instability for one-planet systems could 
come from tidal orbital decay and potential ly direct engulf- 
ment by the rapi dly expanding st ellar envelope. I Villaver fc Livig 
l|2007l . I2OO9I ) and IViUaverl l|201lD treat this effect in detail with 
additional physics s uch as frictional drag, pla net accretion and 
planet evaporation. iMustill fc Villaveil (|2012l ) model individual 
thermal pulses and demonstrate how they affect planetary sta- 
bility. In this study, we only consider planets that are too distant 
to be affected by the stellar envelope expansion and so are not 
affected by tides, accretion or evaporation (see Section |3}. 

1.2.2 Two-Planet Instability 

iDebes fc SigurdssonI (|2002l ) considered adiabatic evolution of 
two-planet systems while exposed to a f Mq star losing half of its 
mass over 1000 planetary orbits. This foundational study consid- 
ered circular and coplanar equal mass planets, with planet/star 
mass ratios ranging from 10~^ to 10~^. They discovered impor- 
tantly that although adiabatic evolution causes both planets to 
move outward and maintain their initial semimajor axis ratio, 
their critical Hill separation changes. 

The rate of change of the separation measured in units of 
Hill's radii is equal to pr'^^^dp/dt. This dependence causes pre- 
viously Hill stable planetary systems to become unstable, and 
incite gravitational scattering which could not occur on the MS. 
Their simulation results suggest that scattering instabilities may 
be more widespread during post-MS evolution than during MS 
evolution. Here, we investigate this claim in significant detail. 
For discussion on the hig h mass-loss non-adiba tic multi-planet 
case recently presented bv lVovatzis et all l|2013l ). please see Sec- 
tion 6.4. 

1.3 Paper Outline 

We begin in Section 2 with a description of the challenges of 
using N-body numerical simulations for gravitational scattering 
amidst mass loss. In Section 3 we determine the regimes where 
engulfment and tides from stellar envelope expansion can be 
neglected for this study. Section 4 presents a general formulation 
of the Hill stability limit and shows how it changes due to stellar 
mass loss. We use the results of Sections 2-4 to motivate the 
setup for our numerical scattering simulations. In Section 5, we 
perform these simulations, and report the results. We discuss 
the consequences in Section 6 and conclude in Section 7. 



2 NUMERICAL CONVERGENCE 

In this section, we highlight the difficulty in achieving accurate 
N-body simulations that model both central star mass loss and 
gravitational scattering amongst multiple massive planets, and 
implement a solution. 




-8 -9 -10 -11 -12 -13 -14 -15 -16 
Logio(Tolerance) 

Figure 1. Difference between interpolating S5-B-outputted mass at 
every Mercury timestep (blue circles) versus interpolating this mass 
within Mercury timesteps (orange squares). "Spliced" indicates the 
latter and "Not Spliced" indicates the former. Shown are the final 
values of the inner planet's semimajor axis for a pair of O.OOIMq- 
planets with initial semimajor axes of 10 AU and 30 AU, and initial 
eccentricities of 0.0 and 0.5, respectively. The parent Solar-metallicity 
star was modeled to lose 6.22Mq of its initial 7.66M0 at the start 
of an asymptotic giant branch phase lasting almost 5 X 10^ yr. The 
"Analytical" line refers to the final semimajor axis of the inner planet 
predicted by adibatic mass loss. The convergence properties of the 
spliced BS integrator for systems with both gravitational scattering 
and mass loss is a significant improvement. 

2.1 Stellar Evolution Code 

We utilize the SSE stellar evolution code (|Hurlev et al.ll2000l l. 
which adopts empirically-derived algebraic formulations in or- 
der to quickly generate a stellar evolutionary track solely from 
a given progenitor mass and metallicity, and stellar model pa- 
rameters such as the Reimers mass-loss coefficient. We use the 
same mass loss pr escriptions as described in Section 7.1 of 
iHurlev et al.l (|200(]| ) with their Reimers mass-loss coefficient de- 
fault value of 0.5. Their choice is observationally motivated 
by the Horizontal Branch m orphology in Galactic globular 
clusters l|lben fc Renzinil Il983|). and lies in the cen ter of the 
range recently considered by Veras fc Wvatj (|2012l '). who dis- 
cuss this choice in light of a n updated version of the Reimers 
law (jSchroder fc Cuntzll2005l ). The SSE code allows us to sam- 
ple many different evolutionary tracks easily, and outputs the 
important parameters, Mi,(t) and -R*(t), where J?j, is the radius 
of the star. 

2.2 Planetary Evolution Code 

We also use the Mercury integration package (|Chamberslll999h , 
which specializes in modeling planetary dynamical evolution. 
In order to accurately model close encounters between plan- 
ets and the parent star - a necessity for this study - we use 
the Bulirsch-Stoer (BS) integrator. This integrator features an 
adaptive timestep, which is determined by a tolerance param- 
eter given at the start of the simulat ion. A tolerance of 10~^ ^ 
is considered to be highly accurate (jjuric fc Tremain3 12OO8I ) . 
Smaller toleranc e values should roughly converge to the same 
result; Fig. 7b of lSmith fc Lissauej (|2009l ) demonstrates that in 
crowded 5-planet systems separated by several Hill radii, toler- 
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ances of 10"^% 10"", 10"^", 10"" and 10"" will yield instabil- 
ity timescales which are all within the same order of magnitude. 
Tolerances below 10"^® generally cannot be achieved because 
in those cases the accuracy requested is greater than machine 
precision. 

2.3 Merging Both Codes 

IVeras et all l|201lf ) found that linearly interpolating SSE stel- 
lar mass output at each Mercury timestep adequately models 
the dynamical evolution of a single planet amidst stellar mass 
loss because the numerical simulations reproduced the analyt- 
ical results. For multi-planet systems, this technique alone is 
inadequate. The interaction between both planets coupled with 
stellar mass loss causes a failure of convergence of orbital pa- 
rameters as the tolerance is decreased. 

In order to improve the accuracy, we have performed an ad- 
ditional interpolation of the SSE stellar mass output in between 
each Mercury timestep at each BS substep. The resulting finer 
gradation makes a crucial difference, as demonstrated by Fig. 
[1] The figure plots the final semimajor axis values for the inner 
10"'^ M0 planet and the outer 10"'^ M© planet in a system with 
initial semimajor axes of 10 AU and 30 AU, and initial eccen- 
tricities of 0.0, 0.5, respectively. All initial orbital angles were set 
to 0°. The simulations were run for the entire evolution of the 
Thermally Pulsing Asymptotic Giant Branch (TPAGB) phase 
of a Z=Z0 = 0.02 (Solar metallicity) star. We chose a progen- 
itor mass of 8M0 to model particularly violent mass loss. Our 
simulations ran during the TPAGB phase only, when Mi, was 
reduced from 7.659M0 to 1.438M0 in about 492,744 yr. 

The plot contains two curves from the simulation output, 
representing final values of the semimajor axis due to SSE- 
outputted mass interpolation at each Mercury timestep (blue 
circles; "non-spliced"), and with an additional interpolation in- 
between timesteps (orange squares; "spliced"). The third, green, 
curve is the analytic prediction for the final semimajor axis of 
the inner planet. This value can be determined because the mass 
loss is adiabatic and the planets are not near a strong mean 
motion commensurability (hence their semimajor axes remain 
secularly unaffected). 

Without the additional interpolation, the results do not ap- 
pear to converge until perhaps at the machine precision limit 
for the BS toleranc^. Further, the extent of the variance in the 
non-spliced curves may fundamentally change the endstate of 
the system if any more close encounters occur. Therefore, we 
use the spliced BS integrator throughout the rest of this work. 
Convergence with the spliced integrator is achieved in this case 
at an accuracy of ~ lO"^'^; we are conservative and adopt the 
value of 10"^'^ for our integrations. 

Another consideration is the ejecta-crossing lag time. Stel- 
lar ejecta will cause the inner planet's orbit to shift before the 
outer planet's orbit. In some cases, this "lag time" between or- 
bital shifts may produce a noticeable change in the dynamics 
that is missed by assuming both planets simultaneously change 
their orbits. The weakness of this assumption is accentuated for 

^ In the one-planet case, when the outer planet is removed from these 
particular simulations, then all three curves are visually indistinguis- 
a ble from one anothe r on this plot. This result reinforces the finding 
of IVeras et al.l l|201lh that splicing within timesteps is generally not 
necessary in one-planet simulations. 



widely-spaced orbits and for systems which are not in the adia- 
batic regime. For the (adiabatic) systems studied here, however, 
this assumption likely produces a negligible effecl[f|, and hence 
is neglected for the remainder of this study. 

2.4 Further Adaptations 

All orbital elements in this work are reported in Jacobi coordi- 
nates. Therefore, as Mercury receives input in astrocentric coor- 
dinates, we performed the conversion. Further, we had to modify 
the default version of Mercury to account for a changing stel- 
lar mass in the output file xv.out so that the conversion from 
Cartesian output to Jacobi elements was performed correctly. 
Consequently, the size of xv . out nearly triples in size. Although 
this increase might be prohibitive for high-resolution studies of 
individual systems, here we are interested primarily in the final 
stability state of each system. Therefore, in our case, outputs 
at a Myr resolution are all that is required. Independent of the 
paucity of outputs. Mercury does record the times of collisions 
to within a timestep. 



3 TREATING THE STELLAR RADIUS 

Additionally, we modified Mercury to incorporate the stellar ra- 
dius evolution profiles from SSE. Over its lifetime, a star's ra- 
dius evolves nonlinearly and nonmonotonically. Because these 
variations are modest and all occur within 0.05 AU during the 
MS, most previous planet scattering studies treated the radius 
as static and/or negligible. However, during the post-MS, the 
radius variations can be violent and extend beyond several AU. 

3.1 Expansion 

Because of the potential for planetary collisions, evaporation 
and/or envelopment due to the expanding stellar envelope, the 
variations in stellar radius must be taken into account dur- 
ing post-MS scattering simulations. An important question is 
whether or not a planet, expanding its orbit due to mass loss, 
can outrun an expanding stellar envelope. The answer is com- 
plicated by the fact that the timescales for and amplitudes of 
mass loss and radial increase are not completely in sync, al- 
though they often are similar. Additionally, a star's radius may 
decrease. WD radii are even more compact than MS radii. 

In order to characterize these variations, we generated Figs. 
OH with SSE data. The figures characterize the maximum stel- 
lar radius for different metallicities and stellar phases, respec- 
tively. In particular. Fig. [2] suggests that the maximum stellar 
radius is largely independent of metallicity, and that roughly the 
number of AU at the maximum radius is equal to the number of 
initial Mq. Figure [3] illustrates that the stellar radius generally 
increases during post-MS phases, although for the IM0 case, 
there is an order-of-magnitude decrease after the RGB stage. 
This decrease becomes progressively smaller as the progenitor 
stellar mass is increased until vanishing at about 3M0. 

^ By using the observed m ass ejecta speed in the post-MS system R 
Sculptoris of 14.3 km/s (iMaercker et al.ll2012t) , one can estimate 
that the ejecta will take 121 days to travel 1 AU. Thus, for an inner 
planet at 10 AU and an outer planet at 12 or 13 AU, the travel time is 
less than a year, a small fraction of the inner planet's orbital period. 
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These trends are model-dependent. Other stellar evolution 
codes (e.g. [Vassiliadis & Wood 1993) demonstrate that as a re- 
sult of thermal pulses on the AGB the expanding stellar radius 
can expand up to 1 AU more than the rule-of-thumb maximum 
radius from the last paragraph. Further, the maximum radius 
at each phase might differ depending on the model used; see, for 
example, iVillaver fc Livig (,2009. ') . 



3.2 Tides 

The maximum stellar radius is just a physical boundary; stellar 
tides can extend beyond the reach of the star. Because modeling 
tides is both beyond the scope of this study and remains the 
subject of debate, we choose initial conditions for our numerical 
simulations where we can safely neglect tidal effects. Planets 
are not necessarily destroyed by tides nor by being engulfed in 
the stellar envelope. The remarkable sub- 10 hour periods of the 
two planets in the hot B subdwarf star KIC 05807616 system 
l|Charpinet et al.1l201ll ) provide strong evidence that planets can 

survive de ep immersion int o the s tellar envelope. 

Both IVillaver fc Livid l|2009l ) and iKunitomo et all (HoTl) 
have analyzed planet engulfment by Red Giant Branch (RGB) 
stars by including a number of physical factors, and use more de- 
tailed stellar evolution models than SSE. Their results indicate 
that, for lower-mass stars, tides can significantly affect plan- 
ets with semimajor axes that are about 2.5 times as high as the 
maximum stellar RGB radius, and that engulfment is sensitively 
dependent on progenitor stellar mass. However, the RGB radius 
does not greatly exceed 1 AU. Although AGB radii are signif- 
icantly larger (Figs. [3 and [3|, iMustiU fc Villaved (|2012l ) found 
that for higher-mass stars, the most distant Jovian planet that 
becomes engulfed is initially at approximately the maximum 
stellar radius; tides can slightly shrink the orbits of planets for 
about an AU beyond this maximum. Hence, we adopt 10 AU as 
the initial orbit of the innermost planet in our simulations. This 
planet will certainly be safe from engulfment in the envelope, 
and will experience no significant tidal decay except possibly 
from the most massive stars. 



3.3 Radius-based Code Modifications 

Our simulation initial conditions are chosen such that stable 
planetary systems are well beyond the influence of tidal effects. 
However, instabilities which arise during the simulations may 
cause planets to approach or collide with the expanding stellar 
envelope. If a collision occurs, the system is flagged as unstable 
and is stopped. In reality, if a star engulfs a planet, the star's 
mass would increase slightly and as a result its radius might 
change as well. Remaining planets in the system would then be 
affected because angular momentum must be conserved. 

Additionally, Mercury's collision detection algorithm 
mce_cent checks to see if the pericenter of a planet's orbit lies 
within the star's radius. If so, a collision is flagged. However, 
if the star's reflex velocity is sufficiently high, then a collision 
might not occur. Therefore, we modified mce_cent to account 
for the stellar refiex velocity. For planetary mass companions, 
however, the refiex velocities are low, and are not likely to factor 
into collision statistics. 
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Figure 2. Maximum stellar radius as a function of progenitor mass 
and motallicity. Colours denote different metallicities; the blue curve 
(Z=0.02) represents Solar metallicity. Curves end when supernovae 
occur. This plot demonstrates a rule of thumb: generally, the number 
of AU at the maximum radius is approximately equal to the number 
of initial Mq . 
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Figure 3. Maximum stellar radius as a function of progenitor mass 
and stellar phase. The y-axis is logarithmic, and the x-axis is monoton- 
ically increasing in time. Symbols denote different progenitor masses, 
and the stellar phases are MS = Main Sequence, HG = Hertzprung 
Gap, RGB = Red Giant Branch, CoHe = Core Helium Burning, 
EAGB = Early Asymptotic Giant Branch, and TPAGB = Thermally 
Pulsing Asymptotic Giant Branch. For Solar metallicity, stellar radii 
generally increase monotonically throughout post-MS phases except 
for progenitor masses approximately equal to IMq. 



4 HILL STABILITY 

We now make some analytical stability estimates to identify the 
systems of interest for our N-body runs. We particularly seek 
systems that are likely to be stable on the MS and subsequently 
destabilised during the giant or WD stages. lDonniso"nl (|2011 ) 
provides a formulation of Hill stability in Jacobi coordinates 
which allows for arbitrarily inclined and eccentric orbits. His 
treatment is fully general with one exception: the expression for 
the system energy is the two-body approximation. This approx- 
imation is necessary in order to obtain an analytically tractable 
(but not strictly correct) solu tion; the intra ctable terms appear, 
for example, in Eq. (2.27) of lVerasI (|2007h . which provides the 
complete expression for the energy of a three-body system in 
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Hill Stability Curves: i„,=0°, Mi-Mi^Mj 
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Hill Stability Curves: ei„=eout=0, Mi=M2=Mj 
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Figure 4. Critical Hill semimajor axis ratios as a function of stellar 
mass for different eccentricities {upper panel), different mutual incli- 
nations {middle panel) and different planetary masses {bottom panel). 
The purple stars on the bottom of each plot, from right to left, rep- 
resent the eventual WD mass [in brackets in the following] for MS 
progenitor masses of 8Mq[1A4.Mq], 7Mq[1.29Mq], 6Mq[1.14Mq], 
5Mq[1.00Mq], 4Mq[0.87Mq], 3M0[O.75Mq], 2Mq10.64Mq] and 
1Mq[O.52M0]. The figure illustrates that the Hill radius is sensitively 
dependent on eout, im, Mi, and M2, but weakly dependent on 
Mir. For Jovian and terrestrial-mass planets, stellar mass loss can 
change the Hill stability limit by at most a few tenths in aout/ain- 



terms of Jacobi orbital elements. See Subsection 16.51 of this pa- 
per for further discussion on this point. 

In the following, the subscripts "1" and "2" refer to 
the inner and outer planets, the subscript "in" refers to the 
star/innermost-planet binary, and the subscript "out" refers to 
the outer planet properties measured with respect to the in- 
ner binary. Let im represent the mutual inclination of the inner 
and outer binar ies. Then the Hill stability curve is given by 
(|Donnisonll201ll '): 



(l + y^) {y^l3^ + 2y^ cos im + I) = 
2Scr {M. + Mi + M2)y^ 



M|(M*-hMi)^(l- 



where 



din 

aoutM*Mi 



(M^MiV^^ /m* + Mi +M2 (l-ef^) 



V M2 J V {M. + Mi)^ (1-e, 
and with (|Donnisonll2"006l ): 



2 1 

out / 



(1) 

(2) 

(3) 



2 (M* + Mi+ M2) 

X (M.Afi + M*M2 (1 + xc.itf + MiM2xli,) (4) 



1 ~t~ Xcvit ^crit / 

2 . ^r 2 



such that X — Xcrit is the unique real solution to the following 
quintic equation: 



(M* + Ml) x'^ + (3M* + 2Mi) x'^ + (3M* + Mi) x^ 
- {3 Ah + Ml) x"" - (3M2 + 2Mi)x = {M2 + Ml) 



(5) 



Care should be taken when choosing the root of the quartic 
equation in Eq. ([1} when solving for y. 

The li terature is replete wi th special-case solutions to Eqs. 
(HI)-® [see lGeorgakarakoJ2008l for a reviewJJ but typically treat 
the masses as static and define a "separation" as a modulated ra- 
tio of the planetary semimajor axes. Equation ([2]) demonstrates 
why. In order to model how the Hill stability curves change with 
mass loss, one need only to evaluate Eqs. (HJ-© at different 
times during a star's evolution. 

In the circular, coplanar, equal-planetary mass case pre- 
sented injOebcs & Sigurdsson (2002), the critical Hill separation 
is shown to vary by a few tenths due to significant mass loss. We 
have undertaken a wider parameter exploration of phase space, 
and have discovered that this result generally holds true for or- 
bits of any eccentricity, inclination and stellar mass as long as 
the planetary masses are at most about one Jupiter-mass each. 
Our results are presented in Fig. O where each panel showcases 
a different parameter dependency. 

The variation in the Hill stability limit due to stellar mass 
loss is often equivalent to several AU for planets which reside 
beyond about 10 AU. Consequently, if planetary packing (e.g. 



We have discovered two typographic errors in t he previous litera- 
ture: T he last quantity of the LHS of Eq. (2.13) in I Veras &: Armitaed 
||2004^ should no t be squared, and the sign in front of A in Eq. (9) of 
iDonnisonI ll201lh is incorrect. 
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iRavmond et al.l l2009l ) produces planets near the Hill stability 
limit during the MS, then post-MS evolution may trigger insta- 
bility. 

Despite these considerations, one must remember that the 
Hill stability criterion is a sujficient but not necessary condition 
for the planets to remain ordered. Hill stable planets may be 
Lagrange unstable, and planets failing to satisfy the Hill stability 
condition may nevertheless be stable. Regarding the latter case, 
one outstanding feature of Fig. 3] is that moderately eccentric or 
inclined orbits yield critical semimajor axes ratios that are high 
- much higher, for example, than the semimajor axis ratios of 
any adjacent pair of planets in the Solar System. 

Large critical semimajor axis ratios may strongly influ- 
ence the location where two planets become Lagrange stable. 
IVeras fc Armitag^ (|2004l ) show that as the mutual inclination 
of the circular orbits of two equal mass planets increases, the 
critical Hill stability limit becomes a progressively worse indi- 
cator of the separations at which planets may actually become 
Lagrange stable. Their Fig. 5 illustrates that for im = 35°, in- 
stability occurs effectively for values under half of the critical 
separation. However, their numerical simulations were run for 
just 2 Myr, almost certainly missing instances of longer-term 
instability. 

Therefore, determining how mass loss affects the stability 
prospects of the orbits of two planets is perhaps more complex 
than just considering the analytic effect on the Hill stability 
limit. Hence, we now turn to N-body simulations. 



5 N-BODY SIMULATIONS 

Ideally, we could sample the entire two-planet/single-star phase 
space with numerical simulations. Realistically, we must take 
measures to restrict our studies to computationally feasible and 
insightful simulations. To better understand how to restrict the 
phase space, we consider typical MS lifetimes, the mass lost at 
each evolutionary phase, and the planetary period enhancement 
during the WD phase, in Subsections 5.1-5.3. We motivate and 
state our initial conditions in Subsection 5.4. Subsections 5.5 
and 5.6 present the simulation results. 
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Figure 5. The MS age of stars with Solar metallicity (blue, top 
curve) and with a very low (Z=0.0001) metallicity (red, bottom 
curve). Higher mass progenitors significantly reduce the CPU time 
needed to integrate planets over a star's entire MS lifetime. 
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Figure 6. The number of orbits taken by an isolated planet around 
a star throughout its MS phase. In each pair of equivalently-coloured 
curves, the top curve is for Z= 0.02 and the bottom curves is for 
Z= 0.0001. The number of orbits decrease with higher stellar mass 
because the decreased MS timescale dominates the shortened orbital 
period. These values are important both physically - to determine 
instability - and computationally - to assess the feasibility of inte- 
grating over the entire main sequence phase. 



5.1 Main Sequence Timescales 

MS ages are given in Fig. [S] These ages may vary by Gyrs 
depending on stellar metallicity, and are at least 1 Gyr long 
for any progenitor mass less than 2Mq. With a MS lifetime 
of over 10 Gyr, a Solar-mass, Solar-metallicity star with orbit- 
ing planets is particularly prohibitive to integrate. This long 
timescale helps explain the uncertainties in long-term evolution 
of the Sol ar System planets (|Kholshevnikov fc Kuznetsovll2007l : 
iLaskar et al. 2011). Figure [6] provides estimates for the number 
of planetary orbits that would be achieved during the MS as 
a function of stellar progenitor mass, for a variety of planetary 
semimajor axes. The curves result from the competition between 
the decreased orbital timescale with increased stellar mass, and 
the decreased MS lifetime with increased stellar mass; the lat- 
ter overwhelmingly wins. Note the number of orbits tails off 
significantly as M*(0) is increased, for all semimajor axes. A 
one order-of-magnitude change in semimajor axis corresponds 
to 1.5-order-of-magnitude change in the number of MS orbits. 



5.2 Post-Main Sequence Phase Properties 

The evolution timescales of the intermediate post-MS, pre- WD 
phases are short compared to the MS timescale. Figure [7] pro- 
vides timescales for each stellar phase, and relates the phase 
to the percentage of the star's original mass lost, for Solar- 
metallicity stars. The plot demonstrates that except for the IMq 
case, most mass is lost during the TPAGB and negligible per- 
centages of mass are lost in the other phases. 

5.3 White Dwarf Period Enhancement 

After the star has become a WD, the star stops losing mass 
and gradually cools down. Compared to its MS mass, the WD 
mass is greatly reduced, and cannot exceed the Chandrasekhar 
Limit of « IAMq. The result is that the planetary period may 
be drastically increased. Assume that the planet's evolution is 
entirely adiabatic. A reduction of the star's mass by a factor 



© 2013 RAS, MNRAS 000.[Tlf24l 



8 



Veras, Mustill, Bonsor & Wyatt 




0.01 0.1 1 10 100 1000 
Phase Duration / Myr 

Figure 7. Correlating mass loss fractions, phase durations and pro- 
genitor masses for stars of Solar metallicity and a Reimers mass loss 
coefficient of 0.5. Each curve contains 8 symbols representing stel- 
lar progenitor masses of 8Mq, TMq, &Mq, 5M0, AMq, 3Mq, 2Mq 
and IMq, ordered monotonically. Green squares represent the Ther- 
mally Pulsing Asymptotic Giant Branch (TPAGB), brown upward- 
pointing triangles the Early Asymptotic Giant Branch (EAGB), 
gray downward-pointing triangles the Core Helium Burning phase 
(CoreHe), red diamonds the Red Giant Branch (RGB), yellow circles 
the Hertzprung Gap (HG), and the blue open squares the MS. Most 
MS mass loss is too small for this plot. Except for the IMq case, the 
most important mass loss phases are the TPAGB, which all last on 
the order of 1 Myr. 



of k will cause the planet's semimajor axis to be increased by 
a factor of k. Hence, the planetary period around the WD is 
k^ times the period around the MS star. Figure [8] plots this 
enhancement factor as a function of progenitor stellar mass, and 
demonstrates both that WD planets perform fewer orbits than 
MS counterparts in the same amount of time (with possible 
implications for scattering) and that WD numerical simulations 
may proceed much more quickly than MS simulations. 




12345678 
Progenitor Stellar Mass /Mq 



Figure 8. The planetary period enhancement factor due to post-MS 
evolution, supposing that the planet has evolved entirely adiabatically 
and is isolated from any other perturbations. This enhancement factor 
is independent of semimajor axis. The top curve is for Z= 0.02 and 
the bottom curve is for Z= 0.0001. 

5.4 Initial Conditions 

5.4.1 Fiducial Choices 

The above considerations lead us to choose an integration dura- 
tion of 5 Gyr and Solar- metallicity progenitor masses of 8M0, 
7Mq, 6M0, SMq, 4M0 and 2,Mq. This mass range extends 
down to th e upper mass-end o f the observed range of exoplanet 
host stars (jSato et al.l l2012"bl ). This combination allows us to 
sample an ensemble of multi-planet systems over every phase 
of evolution, including a substantial sampling (over 4.5 Gyr) of 
evolution in the WD phase. Simulation output occurs at a fre- 
quency of 1 Myr. Performing a statistically significant simulation 
ensemble for IMq stars is well-beyond our available resources; 
for more details on the planetary co nsequences of th e poss ible 
evolutionary tracks of IMq stars, see IVeras fc Wvat^ (|2012l ). 

We adopt one Jupiter mass for the mass of each planet 
(Ml = M2 = Mj), assume the planets are on coplanar orbits 
(im = 0°), and assume they have small but non-negligible MS 
eccentricities (ei(f = 0) = e2(f = 0) = 0.1). These eccentricities 
are low compared to the observed MS values for massive exo- 
planets beyond the tidal circularlization limit, but higher than 
the near-circular orbits predicted by core accretion theory. The 
inner planet is initially set at ai = 10 AU to avoid tides with the 
expanding stellar envelope; over 25 known planets are thought 
to harbor a > 10 AU q Also, this semimajor axis guarantees 
that orbital evolution due to mass loss will be adiabatic unless 
in the presence of an event such as a supernova, which is not 
modelled here. 

We perform 632 simulations per ensemble, where each en- 
semble features a different progenitor mass. In each individual 
simulation, the orbital angles (mean anomalies and longitudes 
of pericenter) of both planets are selected from a uniform ran- 
dom distribution. We adopt 79 values of 02 so that we sample 8 
different sets of orbital angles for each 02 value. The range of a2 
values sampled encompasses both the "chaos limit" and the Hill 
stability limit in order for us to sample many diflferent types of 
dynamical behaviour. 

The chaos limit refers to a maximum semimajor axis ratio 

^ See the Extrasolar Planet Encyclopedia at |http: / / exoplanet. eu/] 
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separation at which mean motion resonances do not necessarily 
overlap. This limit is smaller than the Hill and Lagrange stability 
limits, and represents a fu zzy boundary w ithin which instability 
occurs readily and g uicklv. [ Wisdom! (|l980l ) found the chaos limit 
to be 

for two equa l mass circular- orbit planets. Re cently , 
iQuillen fc Fabeil l|2006l ) and iMustiU fc WvattI (|201^ 
expanded on this result by c onsidering bodies' ec- 
centricity. iMustill fc WvattI (|2012t ) discovered that for 
Cm > 0.21{Mp/Mi,)^^'^ , where Cm is the eccentricity of the least 
massive planet, the chaos limit is (for Mp — Mi = M2): 

We set our minimum value of 02 to be less than the limit 
from the more conservative definition (Eq.[6| for each progenitor 
mass to help ensure that we have a tail of unstable simulations. 
Similarly, we wish to have a tail of stable simulations for large 
separations. Hence, we select a maximum value of 02 that ex- 
ceeds the MS Hill stability limit in each case. 

5.4-2 Additional Simulations 

Motivated by the outcome of our fiducial simulations, we per- 
formed two additional ensembles of simulations (632 simulations 
per ensemble). In both, we adopted a stellar progenitor mass 
of 5Mq. The first case assumed different planetary masses; we 
adopted 1 Earth-mass for each planet (Mi — M2 = Me). In 
the second, we adopted 62 — 0.5 to model a moderately eccen- 
tric outer planet. Doing so yields a much wider Hill stability 
separation (see Fig. [3]) than in the fiducial case. 

5.5 Simulation Results: Fiducial Cases 

5. 5. 1 Overview 

Our goal is to identify instability and when it occurs. We define 
instability as Lagrange instability: if the planets at any point 
are found to achieve a hyperbolic orbit, cross orbits, or collide 
with each other or the star. Hill instability includes just a few 
of these possibilities: planet-planet collisions and crossing or- 
bits. Therefore, Hill stable systems may eventually be Lagrange 
unstable. Those that do will feature ejection of the outer planet 
and/or collision of the inner planet with the star. 

We plot instability time vs. initial semimajor axis ratio in 
Fig.[ni which represents our main result. Dots indicate unstable 
systems. No dot appears for a system that has remained stable 
over the entire 5 Gyr integration. If all 8 simulations for a given 
semimajor axis ratio are stable over 5 Gyr, then we place a 
blue star at 10^" yr in the appropriate horizontal position, even 
though the vertical position of the star has no physical meaning 
and is selected for visual impact. 

The figure includes the six ensembles of simulations with 
different progenitor stellar masses. Each p ost-MS phase c hange 
occurs at different times on each plot. See iHurlev et al.1 l|2000l 'l 
for detailed physical descriptions of each phase. Although the 
horizontal lines are close together in Fig. [9l they are clearly dis- 
tinguishable on the zoomed-in Fig. llUl plots. Different progenitor 
masses also cause differences in the location of the Hill stability 



limit. On each plot is a black vertical dashed line, represent- 
ing the Hill stability limit computed from the star's MS mass. 
The black dotted lines represent the Hill stability limits com- 
puted with the star's WD mass (see the Fig. 2] caption for these 
values). 

5.5.2 Physical Description of Figure\^ 

The Hill Stability limits and the onset of post-MS evolution pro- 
vide boundaries within which one can understand the different 
regions in the Fig. [5]plots: 

1) During the MS, represented by the region under all the 
coloured horizontal lines, dots appear predominately inside of 
the MS Hill stability limit and predominately at times under 
10^ yr. Hence, the limit is useful for identifying short-term in- 
stability. 

2) Some dots appear on the MS but outside of the MS Hill 
stability limit in the plots with progenitor stellar masses of SM© , 
4i\f0 and 3A/0 . All these dots indicate long-term instability (oc- 
curring after ~ 10^ yr). The long-term MS unstable simulations 
with initial separations exceeding the Hill stability limit must 
be (and indeed are) Lagrange unstable: where the outer planet 
is ejected and/or the inner planet collides with the star. As 
the progenitor mass is decreased, the number of unstable MS 
systems beyond the MS Hill stability boundary increases. One 
possible reason is because the longer MS lifetimes (see Fig. [S| 
translate into more orbits for the planets (see Fig. [6]), and hence 
a greater opportunity for instability to occur. Another poten- 
tial reason is that smaller planet-star mass ratios broaden the 
boundaries between Hill and Lagrange instability. 

3) MS instability beyond the MS Hill stability limit appears 
to extend only as far as the WD Hill stability limit. However, 
this must be coincidence - due to our choice of initial conditions 
- as the planetary system has no knowledge of the post-MS mass 
loss that will occur. 

4) The WD Hill stability limits ensure that any WD insta- 
bility occurring beyond this limit must be Lagrange instability. 
Our simulations corroborate this statement. 

5) Each plot in Fig. [5] demonstrates that all sys- 
tems become Lagrange stable for 02/01 > 1.55, just in- 
side the 2:1 commensurability. Reinforcing this estimate are 
blue stars which were excluded from the plot (for vi- 
sual clarity) that extend all the way out in an unbroken 
chain to a2/ai ^ [1.779,1.775,1.770,1.765,1.757,1.747] for 
M©(0) = [8Mq,7M0,6M0,5M0,4Mq,3M0]. This sampled 
range helps to establish the robustness of the Lagrange stabil- 
ity boundary for our chosen 5 Gyr integration duration. This 
boundary lies at a distance corresponding to approximately 
[178%, 176%, 176%, 172%, 167%, 163%] of the MS HiU stability 
limit and [138%, 134%, 133%, 130%, 126%, 123%] of the WD Hill 
stability limit. 

In all our cases, planets with initial separations that exceed 
the 2:1 commensurability are stable throughout the 5 Gyr inte- 
gratior(f|. However, the Lagrange stability limit may be higher 
than the values reported here for progenitor masses lower than 
those sampled here. Especially for IMq stars, the outcome is 
uncertain, given the long MS lifetime and strong mass loss over 

® The strong, first-order 2:1 mean motion commensurability perhaps 
plays a role in establishing the Lagrange stability boundary (see, e.g., 
iBarnes fc Greenberell2007l ) for our fiducial cases. 
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both the RGB and AGB (see Fig. [71). The distribution of unsta- 
ble systems in Fig. shows a wide variation of instability times, 
and suggests that longer simulations, perhaps out to the age of 
the Universe, could yield additional instability. 



atively infrequently: [4.0%, 5.4%, 6.9%, 6.7%, 8.7%, 4.6%] of aU 
unstable systems for each stellar progenitor mass. WD instabil- 
ity is not limited to systems with initial separations beyond the 
Hill stability limit. 



5.5.3 Description of Figures\WST^ 

Post-MS pre- WD phase changes can prompt instability, which 
can be seen more clearly in Fig. [TOl The phases which cause 
the greatest mass loss (see Fig. [7| also tend to be the phases 
which are most likely to trigger instability. This tendency is in- 
dicated by the number of dots between the TPAGB and WD 
lines versus the number of dots appearing below the TPAGB 
line. Note however, that the relatively long length of the core 
helium burning phase helps to increase the number of unstable 
systems during that phase. The abundance of dots just above 
the WD line indicates that the TPAGB can unsettle stable sys- 
tems enough to cause slightly delayed instability. In Fig. IIUI the 
WD Ifill stability limit appears to provide an effective boundary 
beyond which post-MS pre- WD instability does not occur. How- 
ever, this apparent boundary again must be coincidence because 
the systems are unaware of forthcoming post-MS mass loss. 

In support of the above claims, we quantify the types of 
instability in Fig.[TT]for each progenitor mass. The figure shows 
six normalized bar plots. The blue, or topmost, bars represent 
the fraction of systems (out of 8) that feature a collision be- 
tween both planets; orange bars represent a collision with the 
central star; yellow bars represent any other type of instability 
(including ejection or periodic instances of a planet attaining 
a hyperbolic orbit); green, or bottommost, bars, indicate sys- 
tems which remained stable for 5 Gyr. Therefore, Lagrange un- 
stable systems are represented by orange and yellow bars. The 
predominance of the yellow bars towards the right sides of the 
plots indicate that any type of collisions become less likely as the 
initial planet separations are increased. Any blue bars beyond 
the MS Hill stability limit indicate planet-planet collisions dur- 
ing the post-MS, importantly demonstrating that after leaving 
the MS, planets are not restricted to (but still predominately 
experience) Lagrange instability. There are no blue bars that 
exceed the WD Hill stability limit, as expected. The height of 
the green bars around the 3:2 commensurability demonstrates 
how it helps stabilize the simulated systems. 

Also of interest is the evolutionary phase at which insta- 
bility occurs. The collision of a planet with a star has been 
proposed to explain both the existence of extreme horizontal 
branch stars without stellar binary companions - as the enve- 
lope of the progenitor gianl^ could be removed by the planet 
l|Charpinet et al.l |20H| : [Bear fc Sokei) [20]J) - and the enrich- 
ment of L ithium seen in a few per cent of stars at all parts of 
the RGB l|Lebzelter et al.[[2012l l. Of particular interest here is 
the planet candidate proposed orbiting the Lithium-rich giant 
BD-l-48 740, whose eccentric orbit suggests a p ast strong scat- 
terin g interaction such as we are considering ([Adamow et al.l 
[2012h . 

Although the fraction of post-MS pre- WD instability can be 
deduced from Figs. l9l and 1101 we have created a separate figure. 
Fig. 1121 which better visualizes the result. Figure [T^] displays 
the fraction of systems which are stable (green bars), and those 
which incur instability on the MS phase (purple bars), on the 
WD phase (orange bars) and in between the MS and WD phases 
(gray bars). Instability during a giant branch phase occurs rel- 



5.5.4 Potential Resonance Behaviour 

Mean motion commensurabilities, shown on the upper x-axes 
of Figs. I9I12I appear to play an important role in affecting sta- 
bility. At these locations, stability may be either enhanced or 
additionally disrupted. A few suggestive instances of these com- 
mensurabilities making a contribution include the 5:4 location 
in the Mt(0) = 8Mq simulations and the 5:3 location in the 
M*(0) = 3Mq simulations. The 3:2 commensurability demon- 
strably provides a protection mechanism for planetary systems 
in each plot. The 4:3 commensurability seems to yield all pos- 
sible outcomes for M*(0) ^ 5Mq, but no stable systems for 
M*(0) < 5Mq. 

Hence, exploring the resonant character of these systems is 
of potential interest. However, our 5 Gyr simulations are not 
well-suited to determine if a given system is locked into mean 
motion resonance because our output frequency of 1 Myr is over 
4 orders of magnitude greater than a typical orbital period. The 
sudden and drastic changes in eccentricity and inclination which 
can arise from purely 3-body interactions (not including any 
type of di ssipation nor ext ernal forces) may act well within 1 
Myr (e.g. [Naoz et al.[[20lil ). and hence disrupt and/or create 
resonances. Additionally, resonance behaviour may manifest it- 
self only periodically due to repeated separatrix crossings, which 
yield different intervals o f libration and circulation of one or 
more resonant angles (e.g. [Farmer fc Goldreich|[2006[ ). Recently, 
[Ketchum et al.[ l|2012l ') has classified this behaviour as "nodding" 
and analytically characterized it. 

Despite these caveats, we have considered the evolution 
from selected resonant angles from our output. Identifying 
resonant systems requires defining a libration centre, maxi- 
mum libration amp l itude about this centre, and a duration. 
[Veras fc FordI l|2009l . [2OI0I ) demonstrated that fully character- 
izing potential resonant behaviour for two massive planets may 
require the sampling of several libration centres, as well as com- 
puting a mean absolute deviation or root mean squared devi- 
ation about each centre, for each resonant angle. Here we do 
not pursue any such analysis, but instead simply point out that 
the majority of stable fiducial systems over 5 Gyr do not ap- 
pear to exhibit resonant libration. These systems include planets 
with Hill unstable separations. Conversely, a smaller number of 
systems do appear to exhibit resonant libration, typically close 
to the strong first-order mean-motion resonant commensurabil- 
ities. This result is expected, as forming resonances of Jovian- 
mass planets from u niformly-sa mpled orbital angles should be 
infrequent (see, e.g.. [Ver3[2007[ ). 

A perhaps better measure of the effect of mean motion com- 
mensurabilities is by considering the geometric mean of insta- 
bility times for each semimajor axis ratio set of 8 simulations 
that produces at least 1 unstable system. We plot these mean 
times in Fig. 1131 where the 4:3 commensurability is shown to 
have a clear effect. 
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Figure 9. Instability times vs. initial semimajor axes ratios for full-lifetime 2-planet simulations. Dots indicate individual unstable systems, and 
blue stars indicate all 8 systems at that separation ratio were stable over 5 Gyr. Blue stars not shown extend out to at least 02/01 = 1.747 in an 
unbroken string in each plot. Each coloured horizontal line represents a stellar phase change, and is at a different position on each plot. The two 
vertical lines represent the Hill stability limit for the MS (dashed) and WD phase (dotted) for each progenitor mass. Any unstable systems on 
the MS beyond the MS Hill stability limit (such as for the 5Mq, AMq and 3Mq cases) or during the WD phase beyond the WD Hill stability 
limit are Lagrange unstable. The plot demonstrates that instability during the WD phase can be achieved at separations that well-exceed both 
the MS and WD Hill stability Umits. 
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Figure 10. Zoomed-in versions of Fig. |9]to sliow detail at times of stellar phase changes. In ascending vertical order, the phases are "Hertzprung 
Gap" = blue; "Red Giant Branch" = red; "Core Helium Burning" = green; "Early AGB" = orange; "TPAGB" = purple; "WD" = gray. The 
longest pre- WD post-MS phase is the Red Giant Branch; the most violent phase (with the greatest mass loss, and causing the greatest amount 
of instability) is the TPAGB. The WD Hill Stability limit acts as an effective empirical boundary for pre- WD post-MS instability. 
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Figure 11. Types of instability. Blue indicates the fraction of systems (out of 8 per bar) that went unstable because of a collision between the two 
planets. Orange represents instability due to a collision with the central star. Yellow indicates any other type of instability, which predominantly 
includes ejection. Green indicates no instability. Hence, Hill unstable systems are included in the blue bars only, and Lagrange unstable systems 
are included in the yellow and orange bars only. The black dashed and dotted lines are the MS and WD Hill stable boundaries, as in Fig. |9] 
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Figure 12. Phases of instability. Green indicates the fraction of systems (out of 8 per bar) that were stable over 5 Gyr, as in Fig. 1111 Purple 
indicates that instability occurred on the MS. Gray indicates that instability occurred between the MS and WD phases. Orange indicates that 
instability occured during the WD phase. The plot illustrates that instability during a giant branch phase does occur, but infrequently. Also, 
WD instability can occur for Hill unstable planets which survive until the WD phase. The yellow dashed and dotted lines are the MS and WD 
Hill stable boundaries, as in Figs. I9I11I 
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Figure 13. The geometric mean of instability times when at least 
one system out of 8 per semimajor axis ratio goes unstable in each 
fiducial ensemble of simulations. The plot suggests that planets near 
strong mean-motion commensurabilities tend to survive for longer 
times before going unstable. 



5.5.5 Survivor Orbit Properties 

Of potential interest to WD pollution investigators and WD 
planet hunters is the properties of planets undergoing Lagrange 
instability during the WD phase. In some cases, the inner 
planet simply collides directly with the WD, creating a di- 
rect pollutant. However, our simulations suggest that these oc- 
curences are rare, occuring [1,0,2,2,0,0] times for Mq(0) — 
[8Mq , 7Mq , 6M0 , 5Mq , 4M0 ,3Mq]. 

In other cases, the outer planet is ejected and the inner 
planet survives on a bounded orbit. A bound, eccentric inner 
planet may induce pollution by scattering comets or asteroids 
close enough to the WD to be tidally disrupted and ingested 
by the WD. To gain insight into how an inner planet survives 
on an eccentric bounded orbit, we use conservation of energy 
and angular momentum. Although energy is not conserved in a 
system with mass loss - and certainly not in our integrations - 
after the parent star has become a WD, mass loss ceases and 
then energy is conserved for the future. Thus, we can compare 
the states at the beginning of the WD and at the moment of 
ejection (= tins). Further, because the mass loss is adiabatic, 
we can relate the semimajor axes of the planets on the MS and 
the WD phases through knowledge of how much mass is lost. 
Therefore, we find that the semimajor axis of the bound planet 
should be at most: 



ai{ti„ 



M^{0)ai{0)a2{0)Mi 
M^twD) [a2(0)Mi -f- ai(0)M2] 



(8) 



where M*(tivD) is the mass of the WD. 

Angular momentum is conserved throughout a planetary 
system's life, even under the effects of mass loss. Thus, in prin- 
ciple, one can use conservation of angular momentum to deter- 
mine the value of ei(fins). However, doing so requires knowledge 
of the hyperbolic values of 02 (tins) and 62 (tins). These values are 
set by the ejection velocity, which is determined by the strength 



of the instability in each case. Our numerical simulations show 
that ei(tins) varies considerably. 

We plot the semimajor axes (blue squares) and periastra 
(orange dots) of the surviving planets in systems featuring ejec- 
tions in the WD phase only, in Fig. [JH Superimposed on the 
plots through aqua lines are the analytically predicted maximum 
values of ai(tins) through Eq. (|8}. Also plotted as a solid black 
horizontal line is the maximum stellar radius (see Figs. [2] and 
|3} attained during the star's evolution. The presence of orange 
dots below the black line suggest the presence of a population 
of highly eccentric planets orbiting WDs whose present pericen- 
tres take them inside the maximum AGB radius. These planets 
cannot have been formed in situ or anywhere near their WD 
locations because otherwise they would have been destroyed or 
suffered radical orbital alterations on the AGB. The fractions of 
orange dots below the black line for each progenitor mass are 
8.4% (SAfe), 10.9% (7M0), 8.4% (6M0), 5.1% (5M0), 1.7% 
(4Mq) and 0% (SMq). 



5.6 Simulation Results: Additional Cases 

5.6.1 Terrestrial Planet Masses 

Now we consider variations on the fiducial case. First, we set 
the planet masses such that Mi = M2 = Mq. The results of 
those 632 simulations are summarized in Fig. 1151 which include 
plots for instability time versus initial semimajor axis ratio, the 
types and phases of instability, and the geometric mean of insta- 
bility time. The smaller planetary masses here shrink the chaos 
limit and the Hill stability limit - as well as the difference be- 
tween the MS and WD Hill stability limits (see bottom panel of 
Fig. 13} - and introduces a different set of potentially important 
commensurabilities, displayed on the top x-axis of all plots. 

Earth-mass planets fail to go unstable beyond the MS Hill 
stability limit in all but a few cases; the blue stars continue in 
an unbroken chain out to 02/(11 ~ 1.510. Therefore, Hill sta- 
bility and Lagrange stability appear to have almost identical 
boundaries in this case. 

No system features a collision of a terrestrial planet with 
the WD. This result should not imply that this type of collision 
cannot occur, but rather that giant planet collisions are likely 
to be much more frequent. 

Almost all instability involves collision between both plan- 
ets. Interestingly, one planet-planet collision occurs just beyond 
the WD Hill stability boundary, meaning that the real bound- 
ary must differ from the line on the plot. This situation arises 
because planets arrive on the WD phase with slightly differ- 
ent osculating orbital parameters than they harboured at the 
start of the MS (primarily due to their mutual interactions, and 
slightly due to post-MS mass loss). Hence, the real WD Hill 
stability limit for each individual system is different, and dif- 
fers from the line shown in the figure that was computed for 
ei = 62 = 0.1 exactly. 

Although only one system undergoes instability in the post- 
MS pre- WD phases, instability during the WD phase is common, 
and is in fact greatest for separations close to the tightly-packed 
11:10 commensurability (corresponding to 02/01 « 1.066). 
Other instances of WD instability occur around the 7:6 and 
5:4 commensurabilities, and to either side of the 4:3 commensu- 
rability, on which all 8 systems are stable. The influence meted 
by these first-order commensurabilities therefore appears to be 
extensive. However, the bottom plot of the figure might suffer 
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Figure 14. Orbits of inner planets that survive the ejection of the outer planet during the WD phase. Blue squares represent semimajor axes, 
and orange dots represent periastra, all from the simulation outputs. The aqua curve is the analytical estimate for where the maximum semimajor 
axis should be (Eq. (S]!. The black curve indicates the maximum stellar radius achieved during its evolution. Dots below this curve indicate the 
existence of WD planets orbiting inside the maximum AGB radius even though the planets were formed elsewhere. 
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bias due to small number statistics. An investigation of indi- 
vidual systems indicates that resonant behaviour appears to be 
more common for these terrestrial-mass planet systems than for 
the fiducial Jovian-mass systems, perhaps demonstrating the 
importance of the planet/star mass ratio in creating resonance 
behaviour from random initial conditions. 



5.6.2 Different Planetary Eccentricities 

Next, we consider systems with an initially moderately eccentric 
outer planet orbit (e2 = 0.5) and the same fiducial inner planet 
eccentricity (ei = 0.1) in Fig. 1161 In this case, the analytic MS 
Hill stability limit is much larger (02/01 ~ 2.73) than in the 
fiducial case, which is expected (see Fig. 3]). 

For these systems, instability beyond the MS and WD Hill 
limits is extensive. In fact, the Lagrange unstable region extends 
out nearly 7 AU from the location of the MS Hill stability bound- 
ary for tti = 10 AU. However, the Lagrange stability boundary 
lies at a distance that is 143% and 135% of the MS and WD Hill 
stability limits; the former value is a much smaller factor than 
in the fiducial cases. 

Collisions with the star are infrequent and ejections are 
common, perhaps because the inner planet has the much lower 
eccentricity. In fact, inner planets which survive outer planet 
ejections largely fail to attain high enough eccentricities to in- 
trude within the maximum stellar radius achieved during the 
host star's evolution: only 0.85% of the orange dots are below 
the black line in the bottom-right plot. In that same plot, the 
two outliers likely represent systems that experienced pre- WD 
instability that was missed by our low output frequency. Insta- 
bility during the giant phases is uncommon, and is restricted to 
the region around the 4:1 commensurability. 



6 DISCUSSION 

6.1 Consequences for White Dwarf Pollution 

6.1.1 Background 

One potential consequence of dynamical instabilities induced 
due to stellar mass loss is WD pollution. WDs are surrounded by 
thin atmospheres of either hydrogen or helium. Heavier elements 
sink rapidly in such thin atmospheres and it is therefore puzzling 
that such a high fraction of WDs have evidenc e for metal pollu- 
tants in their spectra (25% of single DA WDs. lZuckerman et al.l 
I2OO3I ). Such metal pollution has been associated with excess 
emission in the infrared c onsistent with a c lose-in dust disc in 
more than a doz en cases ( Farihi et al.l 200S ) and gas discs in a 
handful of cases (|Gansicke et al.ll2006l . l20OT . |2008| ). It has been 
suggested that these observations could be explained by plan- 
etary material, motivated by the similarity of the composition 
of the accreted materia l with planetary material l|Klein et al.l 
I2OIOI : iGirven et aLll2012l ) . 

In order for a WD to be polluted by material from an outer 
planetary system, comets, asteroids or planets must be scattered 
at least close enough to the star (at about one Solar radius) so 
that they are tidally disrupted. Changes to the dynamics of 
the planetary system following stellar mass loss has been sug- 
gested as a potential cause of increased n umbers of planetary 
bodies sc attered o nto star-grazing orbits (IDebes fc SigurdssonI 
l2002l : [jur a 2008; B onsor et al]|201ll : lDebes et al.ll2012l '). Even in 
planetary systems where the planets remain on stable orbits. 



iBonsor et al] (|201ll ) and iDebes et al] l|2012l ') show that suffi- 
ciently many asteroids or comets can be scattered onto star- 
grazing orbits to explain some of the observations of polluted 
WDs. Such mechanisms, however, struggle to ex plain observa- 
tions of high accretion ra tes in old polluted WDs (|Koester et al.l 
l201ll : lGirven et al.|[20ll ). 

Instabilities in planetary systems could provide a potential 
explanation for pollution in these, and other, WDs. Depending 
on the exact nature of the instability and structure of the indi- 
vidual planetary system, in many cases the number of planetary 
bodies scattered onto star-grazing orbits increases significantly 
following an instability. This means that the planetary systems 
presented in this work in which instabilities occurred during the 
WD phase have the potential to produce polluted WDs. 



6.1.2 Our Simulations 

We qualify the following discussion by reminding the reader that 
we have not probed progenitor masses between IMq and 2Mq, 
where the true WD population appears to peak. The reason for 
not considering this range is due to computational limitations, 
as tracking planetary orbits over such long MS lifetimes is pro- 
hibitive. 

Figure [5] showcases the different types of instability which 
can occur during the WD phase. First, many systems with giant 
planets that were Hill stable on the MS become Hill unstable 
due to mass loss preceding the WD phase. Second, other sys- 
tems that were technically Hill unstable on the MS, but were 
protected against instability by commensurabilities, become un- 
stable during the WD phase. Third, some of the systems that 
are Hill stable during the WD phase are actually Lagrange un- 
stable, and experience instability at a late time. 

We can now relate the WD instability to the cooling age of 
the WD, which is the time lapsed since the star became a WD. 
The cooling age of real (not simulated) WDs are readily deter- 
mined from observations, and hence provide an opportunity for 
comparison with and motivation for numerical simulations. Fig- 
ure [T7| displays the number of planetary systems that become 
unstable during the WD phase as a function of cooling age; each 
plot corresponds to a different progenitor mass. Any instance of 
instability could represent a trigger for a potentially polluting 
event. An ejection especially suggests that the inner planetary 
system will be significantly perturbed, and this perturbation can 
throw planetesimals or terrestrial planets that were originally in 
the inner system onto the star. 

The figure contains several notable features: 1) planet- 
planet collisions tend to occur at short cooling ages (~ 10'' — lO'' 
yr), 2) escape and stellar collisions becomes prevalent only after 
~ 10^ yr, but tails off after a few Gyr, 3) the cooling age of this 
type of instability increases as the progenitor mass increases, 
and 4) direct collisions with the WD are rare. The striking delay 
of ~ 10^ yr before escape instability becomes dominant reflects 
the (largely unexplored) dependencies of Lagrange instability on 
number of orbits completed, the initial values of the semimajor 
axes, and the mass ratios. In fact, perhaps the positive corre- 
lation of cooling age instability with progenitor mass is due to 
more orbits being required for smaller mass ratios for Lagrange 
instability to kick in, a feature also apparent on Fig. [§] on the 
MS. Alternatively, the correlation may result from the much 
wider orbits of the giant planets around higher-mass WDs. 

The few giant planets which hit the WD are particularly 
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Figure 15. Here, Mi = M2 = IM^. The plots are of the same format as in Figs. [9] [10] [11] [12] and [HI Instability during the WD phase 
occurs almost exclusively at initial separations inside the MS and WD Hill stability limits. MS Hill unstable systems protected by mean motion 
commensurabilities, many of which feature resonant behaviour, allow for survival during the MS before becoming unstable on the WD phase. 



interesting as pollution sources, and imply that some polluted 
WDs may result from the disruption of giant planets rather than 
comets or asteroids. Planetary material that is accreted onto 
WDs might be composed of many small asteroids or comets 
l|jural I2OO3I : iBonsor et al]|201ll: lOebes et al.,i2012.) . but could 



also be the result of a single large object. Detailed compositional 
analysis of some objects concludes that the accreted material re- 
sembles more closely the bulk Earth in composition than chon - 
dritic material (e.g. IZuckerman et aL I l2007l : iKlein et aLlboidl : 
it is feasible that a disrupted planet would produce high levels 
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Figure 16. Here, 62(0) = 0.5. The plots are of the same format as in 
Lagrange instabiUty can occur readily in moderately eccentric systems. 
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|9][l0l[ni[l2l[l3]and[l4] This figure importantly demonstrates that 



of pollution in a WD. Such disrupted planets could provide the 
explanation for a handful of polluted WDs, in particular old, 
heavily polluted WDs. However, our simulations indicate that 
star-grazing giant planet collisions are too rare to alone account 
for the abundance of observed differentiated material. This re- 
sult agrees with the implied low fraction of planetary collisions 



with WDs from the finding that at most 1-5% of WDs have high 
accret ion rates due to dust disks l|Debes et al.ll20l"ll : lFarihi et al.l 
l2012d ). 

Observed polluted WDs can provide the inspiration for ex- 
tensions to our work, with numerical simulations tailored to par- 
ticular observational samples or campaigns. However, the pur- 
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Figure 17. Instability as a function of the time elapsed since the end of the of the Asymptotic Giant Branch phase, when the WD was born. 
Planet-planet collisions occur quickly (~ 10'* — 10^ yr), whereas instances of escape or stellar collisions typically do not occur until after a few 
tens of Myr. 



pose of this paper is not to make a detailed comparison with the 
observed polluted WD population, which is heavily biased and 
suffers from uncertainty of the distribution of their progenitor 
masses. Our computationally-expensive simulations were set up 
to demonstrate that full-lifetime orbital integrations including 



realistic mass loss are now achievable, and to explore dynamical 
properties of the resulting instability. Nevertheless, we can make 
a crude comparison of the potentially WD-poUuting events in 
our simulations (Fig. I17p and the observed distribution of pol- 
luted WDs (in Fig [18}. Data for the observed distribution of 78 
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WD a ges was taken from lFarihi et all (|2012bh and lGirven et all 
ll2012D . 

Both the observations and simulations show a broad con- 
sensus. Observations confidentl y detect pollution at c ooling ages 
from a few tens of Myr (e.g. iCansicke et al.l |2012| ) to several 
Gyr l|Koester et al.|[201ll ). but cannot detect pollution at earlier 
times when the WD is too hot to rule out a primordial origin for 
metalfl nor at later times when the WDs are too intrinsically 
faint. The simulations begin to demonstrate ejection instability 
only at times greater than tens of Myr due to a true physical 
effect which may be explained by the timescale for Lagrange 
instability to act on the WD phase. This instability appears to 
have largely run its course after a few Gyr have passed, so that 
relatively fewer systems are likely to become unstable beyond 
the 5 Gyr integration time. Observationally, comparing the pol- 
lution frequency at 5 Gyr with that at tens of Myr is difficult 
due to intrinsic biases, despite the s uggestive nat u re of the dis- 
tribution in Fig. 1181 Nevertheless, iBonsor et al.l l|201lh found 
that asteroidal accretion onto WDs from particles scattered due 
to the overlap of mean motion resonance exterior to the planet 
follows a similar trend of being present beyond a few tens of 
Myr after post-MS evo l ution before eventually tailing off after 
a few Gyr. iDebes et all l|2012D later found that asteroidal accre- 
tion onto WDs from both exterior and interior resonances also 
follow this trend. 



6.2 Phase Space Extrapolation 

The results of this study raise several questions. One important 
question is how robust our conclusions are to variations of the 
initial conditions. Because of the computational expense of run- 
ning simulation ensembles for 5 Gyr with the Bulirsch-Stoer al- 
gorithm, we could not perform a wider phase space exploration. 
However, we can guess how the outcomes would vary in other 
circumstances. 

Varying ai would change the number of orbits completed 
over the main sequence; the effect is similar to varying M*. 
Hence, Fig.[9]demonstrates the likely consequence: if the planets 
complete enough orbits, they will be prone to long-term MS 
instability. If the initial planetary separation is beyond the Hill 
stability limit, then this instability must be Lagrange instability. 
Otherwise, the type of instability is unrestricted. 

The consequence of varying planetary eccentricities and 
mutual inclinations is less clear. What is clear is that planets on 
eccentric orbits will feature in post-MS systems. Observations 
suggest that if giant planets are formed at several AU on circular 
orbits, they are unlikely to retain their primordial eccentricities 
on the MS. According to the Extrasolar Data Explorei[f|, as of 
31 Oct, 2012, only 20.4% of exoplanets with Mp > O.lMj and 
a ^ 1 AU have e < 0.1. This percentage shrinks to 10.4% for 
e < 0.0!^. If one were to include only planets in multi-planet 
systems, these percentages become 18.0% and 14.0%, respec- 
tively. Therefore, if the planets we currently observe predomi- 
nately survive until the post-MS, the significant majority will 



These young WDs are able to radiatively levitate metals, masking 
what could be external pollution. 
* at http://exoplanets.org/ 

^ Admittedly, there is a bias towards fitting higher eccentricity values 
than the true values. The disparity worsens with sparser data, which 
is often associated with the most distant radial velocity planets. 



An Observed Polluted WD Sample 



Nuinber 
of WD 
Polluted 
Systems 




Log^Q [White Dwarf Cooling Ag^ /yr 

Figure 18. Cooling ages for a sample of 78 observed WDs from 
Farihi et al. (2012b) and Girven et al. (2012). Although this figure 
may be compared with Fig. 1171 the progenitor masses of these WDs 
are unknown, and the sample suffers from observational bias. 



enter those phases with nonzero eccentricities. Nevertheless, it is 
of interest to determine if MS Lagrange instability can occur for 
planets formed on perfectly circular orbits. Therefore, we have 
performed 48 additional simulations with ei(0) = 62(0) = 0.0 
for M*(0) = 3Mq at different locations beyond the Hill stabil- 
ity limit. We can confirm that a few of those systems become 
Lagrange unstable. A detailed statistical comparison is best left 
to a more comprehensive phase space study. 

We can estimate the dependence on planetary mass by ex- 
trapolating from Figs. |9] and 1151 The trend suggests that as 
the test particle limit is approached, where there is no mutual 
interaction between secondaries, the difference between the Hill 
and Lagrange boundaries might tend to zero faster than the Hill 
boundary itself. In the opposite limit, for two brown dwarfs and 
a more massive evolving primary, we expect the stability bound- 
aries to widen even further than in our fiducial case. In these 
instances, instabilities might be common, and could represent a 
significant catalyst for WD pollution. 

Adding more planets to model real systems such as HR 8799 
introduces a significant complication, but is a viable avenue of 
future study with our numerical code (Mustill & Veras, in prep). 
If additional planets are of similar masses, then we expect them 
to represent destabilizing influences, particularly in the post-MS 
stages. 



6.3 Comparison with Radial Velocity Exosystems 

The majority of known exoplanets are likely to be engulfed dur- 
ing the post-MS phase, as the population of planets beyond 10 
AU is largely unknown. Nevertheless, scaled versions of many 
known systems could represent genuine test cases for our model, 
because of the observed close separations of pairs of exoplanets. 
Further, our findings of MS instability warrant a closer look at 
these systems. 

We have compiled a list of all pairs of planets in the same 
system with a > 1 AU each that were detected by radial ve- 
locity measurements. We obtained the data from the Extrasolar 
Planets Explorer on 6 November, 2012. Using Eqs. HJ-©, and 
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assuming minimum masses and coplanarity, we determined if 
the systems are currently Hill stable. 

We found that 6 pairs of planets are not (HD 181433 c,d; 
HD 204313 b,d; 24 Sex b,c; BD +20 2457 b,c; HD 128311 b,c; 
HD 200964 b,c). Hence, these planets are likely to be protected 
by their proximity to mean motion commens urabilitie s insid e 
small islands of stability, as in, e.g. IWitten mver et al.l (|2012l ). 
This behaviour is reflective of some of our Hill unstable systems 
in Figs . [9l [151 and [161 BD +20 2457 is notable because its mass 
is about 2.8M0, and may be related to our 3Mq simulations. 
However, both planets in the system are likely instead to be 
classified as brown dwarfs (Mi « 22.7Mj and M2 « 13.2Mj). 
The correspondingly high mass ratios in the system, coupled 
with their close separations {aijax ~ 1.39), perhaps intuitively 
suggest that instability should be imminent. Hence, the exis- 
tence of such a system is stark confirmation that Hill unstable 
systems can remain stable for long periods. The system's fu- 
ture prospects for stability on the MS would require detailed 
modelling. 

Three pairs of planets are Hill stable, but have a separation 
ratio which exceeds the Hill stability limit by no more than 30% 
(HD 37124 c,d; 47 Uma b,c; HD 183263 b,c). Based on our simu- 
lations, the close proximity of these systems to the Hill stability 
limit suggests that they might not be Lagrange stable for the 
remainder of their MS lifetimes. Testing this suggestion would 
require a detailed suite of long-term simulations for each system. 
If Lagrange instability scales strongly with planet /star mass ra- 
tio, then the planets in HD 183263 are perhaps in the greatest 
danger: the mass ratios in those systems are about 9.5 times as 
great as any mass ratio that we considered in our simulations. 

Three pairs of Hill stable planets have separations exceeding 
the Hill stability limit by between 45% and 65% (HD 108874 
b,c; HD 159868 c,b; HD 10180 g,h), and two planet pairs (HD 
4732 b, c; mu Ara b, c) have separations that are over twice the 
Hill limit. Although the Lagrange stability boundary is likely 
dependent on several variables, we have performed additional 
simulations as proof of concept to show that in at least one 
case, this boundary can extend out to twice the Hill stability 
boundary. Hence, HD 108874, HD 159868 and HD 10180 are 
not guaranteed to be Lagrange stable without further detailed 
analyses. 

If the planets in any of these observed systems are not copla- 
nar, then they are more likely to be Hill unstable (see the middle 
panel of Fig. 3)). A mutual inclination of just 12.3° would render 
the planets in HD 37124 Hill unstable. At the opposite extreme, 
for the widely separated planets in HD 4732, a mutual inclina- 
tion of 61.0° would be required. 



6.4 Non-Adiabatic Mass Loss 

Adiabaticity in the two-body problem with mass loss is well- 
defined fe.g. IVeras et al.ll201ll ). If we were to assume the two- 
body definition of adiabaticity for each of the planets in each of 
our simulations, then we can claim that at no time did our stable 
planetary systems (with a\ (0) = 10 AU) approach a regime that 
featured non-adiabatic mass loss. However, if 2 planets were 
stably orbiting at separations of hundreds of AU on the MS, 
then the mutual planet-planet interaction coupled with stellar 
mass loss could yield unpredictable evolution during the post- 
MS. 



Analytic Hill Energy Approximation Quality 
I ■ ■ ■ ■ I ■ ■ ■ ■ I ■ ■ ■ ■ I ■ ■ ■ ■ I ■ ■ ■ ■ I 
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Figure 19. The error of using the two-body approximation to model 
the total energy of the system. The black curve and filled dots indi- 
cate the mean energy error; the blue dashed curve and filled squares 
the median energy error. These curves are drawn for fiducial initial 
orbital parameters described in Subsection l5.4.l1 The orange filled di- 
amond and upwards-pointing triangle represent the mean and median 
energy error for the ei = 0.5 simulation. The purple filled downwards- 
pointing triangle and open dotted circle are for the simulations with 
the pair of Earth-mass planets. These errors correspond to semimajor 
axis differences on the order of 10"^ AU for a 10 AU. 



IVovatzis et aP (|2013l l has recently explored non-adiabatic 
mass loss in the three-body problem. By using Lypaunov charac- 
teristic numbers to create dynamical stability maps, they found 
that non-adiabatic mass loss can cause a stably interacting pair 
of planets to shift to a choatic region of phase space. The result- 
ing instability (manifested by escape or collision) may be latent, 
sometimes not appearing for a time that exceeds the duration 
of the mass loss by a factor of tens. In cases where the outer 
planet escapes, distinguishing whether the instability was trig- 
gered by Lagrange instability, non-adiabatic mass loss, or both 
might require detailed follow-up simulations. 



6.5 Sharpness of the Hill Stability Limit 

No violations of the Hill stability limit have occurred in our 
simulations, despite us using the two-body approximation for 
the energy in the analytical formulation (Eqs.[T][Sl). Nevertheless, 
we can estimate the error in semimajor axis as a result. Figure 
[T9I plots the mean and median energy error incurred by using 
the two-body approximation for the energy of the system. In 
all cases, the energy error is between 0.02%-0.06%. This range 
corresponds to semimajor axis differences of several 10"'^ AU 
for planets with a « 10 AU. Our gradient of semimajor axis 
values sampled for our fiducial simulations exceeded this value 
in each case, helping to confirm why the analytic limit was not 
violated and suggesting that the unknown integration error was 
comparably small, if not smaller. 
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7 CONCLUSION 

Architectures of planetary systems during each stellar phase 
may represent historical tracers of formation and presage future 
evolutionary instability and death. We have performed 5 Gyr 
simulations that consistently treat the dynamics of two massive 
planets and every phase of stellar evolution for a wide range of 
progenitor stellar masses {3Mq - 8Mq). These computationally- 
demanding simulations suggest that stable MS systems are in 
danger of future instability. The zone of danger is wide, reach- 
ing out to 163%-178% of the MS HiU stability limit and 123%- 
137% of the WD Hill stability limit for our low eccentricity 
(ei(0) = 62(0) = 0.1) simulations. The consequences for WD 
pollution may be significant: For example, the inner planet can 
be perturbed onto a highly eccentric orbit which takes the planet 
close to the WD or hits the star directly. 
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